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Abstract 

This hybrid method (FE-DVR), introduced by Resigno and McCurdy, Phys. Rev. A 62, 032706 
(2000), uses Lagrange polynomials in each partition, rather than "hat" functions or Gaussian func- 
tions. These polynomials are discrete variable representation functions, and are. orthogonal under 
the Gauss-Lobatto quadrature discretization approximation. Accuracy analyses of this method are 
performed for the case of a one dimensional Schrodinger equation with various types of local and 
nonlocal potentials for scattering boundary conditions. The accuracy is ascertained by comparison 
with a spectral Chebyshev integral equation method, accurate to 1 : 10^^^. For an accuracy of the 
phase shift of 1 : 10~^ The FE-DVR method is found to be 100 times faster than a sixth order 
finite difference method (Numerov), is easy to program, and can routinely achieve an accuracy of 
better than 1 : 10^^ for the numerical examples studied. 
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I. INTRODUCTION 



The solution of differential equations by means of expansions into discrete variable rep- 
resentation (DVR) basis functions has become very popular since it was first introduced in 
the early 1960's A review can be found in the paper by Light and Carrington [2|, and 
generalizations to multidimensional expansions are also under development jsj. 

Previously the main application of the DVR method was for obtaining bound state en- 
ergies and wave functions. For this purpose the wave function is expanded into a set of 
basis functions, whose expansion coefficients are to be determined. The calculations are of 
the Galerkin type, namely, the Hamiltonian applied to the wave function is multiplied on 
the left by each one of the expansion basis functions, and the result is integrated over the full 
range of the domain of the variable, leading to a set of linear equations for the expansion 
coefficients. The integrals to be evaluated are then approximated by discrete sums over the 
values of the integrand evaluated at the support points times certain weight factors such as 
in the Gauss quadrature methods [4|. 

In the case of the solution of scattering problems the finite element method (FEM) 
has also been developed. In this procedure the radial range is divided into partitions, also 
called elements, and the solution of the wave equation in each partition is expanded into 
basis functions such as "hat" functions, Gaussians, or polynomials of a given order, whose 
expansion coefficients are to be determined. The equations for the expansion coefficients 
are obtained through a Galerkin procedure, and in many cases the integrals over the basis 
functions can be done analytically. The continuity of the wave function from one partition 
to the next is achieved by imposing conditions on the expansion coefficients, as is done for 
example in Ref. Ib]. In the more recent DVR methods the basis functions are Lagrange 
polynomials whose zeros occur at the Lobatto points [?], j^, in which case the quadrature 
is denoted as Gauss-Lobatto, and the basis set of functions is denoted as Lagrange-Lobatto. 
This basis set was first suggested by Manolopoulos and Wyatt [9'], and a extensive review is 
given in Ref. 10| . The main computational advantage of using DVR basis functions is that 
the sum mentioned above reduces to only one term, because the product of two different DVR 
functions vanishes at the support points, and only products between the same DVR functions 
remain. Furthermore, within the approximation of the Gauss-Lobatto quadrature rule, the 
basis functions are orthogonal. Hence the procedure leads to a discretized hamiltonian (A^ x 
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A^) matrix, whose eigenfunctions determine the expansion coefficients and the eigenvalues 
determine the bound-state energies. There are several types of errors introduced by this 
method. One is due to the truncation of the expansion of the wave function in terms of basis 
functions at an upper limit A^, another is due to the approximation of the Gauss-Lobatto 
quadrature described above in terms of discrete sums over the support points. A third error 
is the accumulation of machine round-off errors. These errors have been examined for bound 

Hi, [isl' 0, and it is found that the convergence of the energy 



state energy eigenvalues, 



with the number A^ of DVR basis functions is exponential, and the non-orthogonality error 
becomes small as A^ increases. 

Very recently a combination of the FE and DVR methods has been introduced into 
atomic physics by Rescigno and McCurdy [14] for quantum scattering calculations . These 
calculations use the FEM approach but in each partition the basis functions are Lagrange 
polynomials, and the support points are Gauss- Lobatto. This "hybrid" method, denoted as 
FE-DVR, is now extensively used for atomic physics calculations, such as for multi electron 
density distributions in atoms 15|, for photo-ionizin g cr oss sections with fast photon pulses 



161], [17|, and for atom-atom scattering calculations [18 



to name a few. However, in these 



works the accuracy of the results was not studied in detail. The FE-DVR method is also 
used extensively for fluid dynamic calculations since the 1980's fl9] and also in Seismology 



20|], where it is called spectral element method. 



The main purpose of the present study is to investigate the accuracy of the FE-DVR 
method for the scattering conditions, since all the errors described above (the Gauss- Lo- 
batto's integration error, the truncation errors of the expansions, and the accumulation of 
round-off errors) are still present. In our study a method of imposing the continuity of the 
wave function and of the derivative from one partition to the next is explicitly given, and 
the accuracy is obtained by comparing the results of the FE-DVR calculation for partic- 
ular solutions of a one dimensional Schrodinger equation with a spectral jsij Chebyshev 
expansion method 22], S-IEM. The accuracy of the latter is of the order of 1 : 10~^^, as 
is demonstrated in Appendix A. In our present formulation of the FE-DVR the so-called 



bridge functions used in references 



14| - 18| in order to assure the continuity of the wave 



function are not used, but are replaced by another method. 

In section II the FE-DVR method is described, in section III the accuracy is investigated 
by means of numerical examples, and section IV contains the summary and conclusions. 
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Appendix A contains a short review of the S-IEM method, in appendix B an estimate of 
the accumulation of errors is presented, and in Appendix C some accuracy properties of the 
finite difference Numerov (or Milne's) method are presented. 



II. THE FE-DVR METHOD 

The FE-DVR version of the finite element method differs from the conventional FEM 
in that the basis functions for the expansion of the solution iIj{x) in each partition are N 
"discrete variable representation" (DVR) functions, which in the present case are Lagrange 
polynomials ii{x), i = 1,2, ..N of a. given order — 1, 



defined for example in Eq. (25.2.2) of Ref. 23(], and in section 3.3(i) of Ref. 



These 



functions are widely used for interpolation procedures and are described in standard com- 
putational textbooks. This FE and DVR combination was introduced in Ref. [l^, and has 
the advantage that integrals involving these polynomials amount to sums over the functions 
evaluated only at the support points. In the present case the support points are Lobatto 

n 

points Xj and weights Wj, j = 1,2,..N, defined in Eq. (25.4.32) of Ref. [23], in terms of 
which a quadrature over a function f{x) in the interval [—1, +1] is approximated by 

/ /(a;)da;^^/(x,)^,-. (2) 
-^-1 j=i 

If / is a polynomial of degree < 2A^ — 3 then Eq. ([2]) will be exact. This however is not the 
case for the product of two Lagrange polynomials ^{x)ij{x), a polynomial of order 2(A^ — 1). 
In the Gauss-Lobatto quadrature approximation [7|, 8|], given by the right hand side of Eq. 
([2]), these Lagrange polynomials are orthogonal to each other, but they are not rigorously 



orthogonal [l2| because the left hand side of Eq. ([2]) is not equal to the right hand side. If 
the integral limits are different from ±1, such as J^^ f{r)dr then the variable r can be scaled 
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in that we do not use their 



to the variable x. Our method differs from that of Ref. 
"bridge" functions, but rather insure continuity of the solution and its derivative from one 
partition to the next by using only the Lagrange functions. Since the Lobatto points are 
not evenly spaced, expansion ([2]) converges uniformly, which is a general feature of spectral 
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methods 2l|. A further DVR advantage is that the Gauss-Lobato approximation of the 
integral 

ii{x)f (x) ij{x)dx ~ 6ijWjf{xi) (3) 



-1 

is diagonal in i,j and is given by only one term. The convolution 

^1 .1 

ii{x) / K {x,x') £j{x')dx' dx ^ WiWjK{xi,Xj) (4) 
-1 J-i 

is also approximated by one non-diagonal term only, which is a marked advantage for solving 
nonlocal or coupled channel Schrodinger equations. The kinetic energy integral can be 
expressed in the form 



(5) 



^+e]^{r) = V{r)^{r), (6) 



after an integration by parts. In the above the prime denotes d/dx. The integral on the right 
hand side of this equation can be done exactly with the Gauss-Lobatto quadrature rule ([2]), 
since the integrand is a polynomial of order 2N — 4, that is less than the required 2N — 3. 

For the case of a local potential V with angular momentum number L = the equation 
to be solved is 

dr 

and for a nonlocal potential K, the term V{r)ilj{r) is replaced by K{r,r')ip{r')dr'. The 
wave number k is in units of fm~^ and the potential V is in units of fm'"^, where quantities 
in energy units are transformed to inverse length units by multiplication by the well known 
factor 2m /h^. In the scattering case the solutions il){r) are normalized such that for r — )• oo 
they approach 

ip{r) — i- sin(A;r) + tan((5) cos(A;r), (7) 
and with that normalization one finds 

1 

tan(5) = — — / sin(A;r) V{r) ip{r)dr, (8) 



as is well known 



24| 



The FE-DVR procedure is as follows. We divide the radial interval into Nj partitions 
(also called elements in the finite element calculations [5]), and in each partition we expand 
the wave function into Lagrange functions ^i(r), i = 1,2, ..N 

N 

^(J\r) = ^ c[^hi{r), b[^'> < r < (9) 



i=l 



The starting and end points of each partition are denoted as h^^^ and 62^-*, respectively. We 
define the value and the derivative of the wave function at the end point of the previous 
partition as 

^(^-)(6r^^)=cr^\ (10) 

where cj^~^^ is the last coefficient of the expansion (ED oii)^^-^\ and 

^ f/'-'\h't'') = Y.^-'^i '^-''). (11) 

respectively, where £ -(r) = dii{r)/dr. The result ffTOl) follows from the fact that that 
^4(^2) = for i = 1,2,..N — 1, and lN{i>2) = 1- For the first partition we arbitrarily take 
a guessed value of A^'^^ for the non-existing previous partition, and later renormalized the 
whole wave function by comparing it to a known value. That is equivalent to renormalizing 
the value of A^^\ In finite element calculations continuity conditions of the wave function 
from one partition to the next are also imposed. However, the method described below 
applies specifically to the case that the basis functions in each element are of the DVR type, 



rather than general polynomials of a given order 

By performing the Galerkin integrals of the Schrodinger Eq. over the ii in each partition 

J 



^(r + v-F)^(^))= (12) 



[' ei{r){T + V -k'^)-^^-^\r)dr = 0, i = l,2,..iV (13) 

we obtain a homogeneous matrix equation in each partition for the coefficients cl'^\ i = 
1,2,. .N 

MiJ) ^J) = 0, (14) 

where c^*^^ represents the {N x 1) column vector of the coefficients cl'^\ and where the matrix 
elements of M are given by Mij = {ii{T + V — k'^)ij). Here T = —d^/dr"^. The continuity 
conditions are imposed by transforming the homogeneous equation fll4p of dimension into 
an inhomogeneous equation of dimension N — 2 whose driving terms are composed of the 
function ip and dip/dr evaluated at the end of the previous partition. These conditions are 
given by 

<=</'= cs;-" (15) 



where use has been made of ii{bi) = for i = 2, ..A^, and = 1, and 



dr 



1=1 



These two conditions can be written in the matrix form 



(16) 



Fiia + F12/3 = 7 



(17) 



where 



where 



where 



and where 



Fu 



1 



(J) 



'12 



'1 -2/ 



a 
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f c 



{J) 



Cl 

C2, 



C4 



Cat 



(J-1) 



With that notation Eq. ( IT^ can be written in the form 



(^) 



(18) 



(19) 



(20) 



(21) 



Mn M12 \ /a 
M21 M22/ \B 



0, 



(22) 



where the matrix M^"^) has been decomposed into four submatrices Mn, M12, M21, and M22, 
which are of dimension 2 x 2, 2 x (A^ — 2), (A^ — 2) x 2, and (A^ — 2) x (A^ — 2), respectively. 
The column vector a can be eliminated in terms of /3 and 7 by using Eq. fll7p . 



« = Ff/(-Fi2/3 + 7) 



(23) 



and the result when introduced into Eq. (122]) leads to an inhomogeneous equation for /3 



:-M2iFri^Fi2 + M22)/3 = -M2iFfii7. 



(24) 



Once the vector /3 is found from Eq. (^^, then the components of the vector a can be found 
from Eq. fl23l) . and the calculation can proceed to the next partition. 
If one expresses the inverse of Fn analytically 

^■11'= I '^°). (25) 

^2 ^2 / 



then one finds 



,(^-1) 



^"7=1 J-n (26) 



^2 '^^ ~^ ^2 



and 




= L, , „ . (27) 



Inserting f l25|) into (!23|) one finds that c\^^ = Cn ^lut c^^ is a function of cj^ A^-^ ^\ 
and the vector /3. 

III. ACCURACY 

We have tested the accuracy for cases with angular momentum L = for two local 
potentials Vm and Vws, shown in Fig. [H and for a nonlocal potential K{r, r') of the Perey- 



Buck type 26| . Potential Vm is of a Morse type with a repulsive core near the origin, given 
by 

VM{r) = 6exp(-0.3 r + 1.2) x [exp(-0.3 r + 1.2) - 2] . (28) 
and V\Ys is a short-ranged simple Woods-Saxon potential given by 

Vws{r) = -3.36/{l + exp[(r - 3.5)/0.6]}. (29) 

The coefficients 6 and —3.36 are in units of fm'"^, the distances r are in units of fm, and 
all other factors are such that the arguments of the exponents are dimensionless. These 
potentials are shown in Fig. [H and the respective wave functions are shown in Fig. O 
The choice of these potentials is motivated by the difference in the degree of computational 
difficulty that they offer in the solution of the Schrodinger equation. Potential Vws has 
no repulsive core near the origin and is of short range. Hence the corresponding wave 
function does not have large derivatives near the origin, and needs not to be calculated out 
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FIG. 1: (Color on line) Potentials Morse Vm and Woods Saxon Vws as a function of radial 
distance r. These potentials are given by Eqs. (128]) and (f29l) . respectively. 



to distances larger than 20 /m, where the potential is already negligible, of the order or 
10^^^. By contrast, neither of these two features apply for the case of Vm- In order to obtain 
an accuracy of 1 : 10~^^ the wave function has to be calculated out to 100 fm, as is indeed 
done in the calculation of the bench mark S-IEM solution, and the repulsive core near the 
origin is more difficult to treat. The nonlocal potential K is described in Eq. (3) of Ref. 25| 
together with the Appendix of Ref. 3] • The accuracy of the corresponding wave function 
obtained with the S-IEM method for this nonlocal potential is illustrated in Fig. 7 of Ref. 



25| . For the nonlocal case only one partition is used in the FE-DVR method, that extends 



from r = to -Rmax? but in view of Eq. (jlj), the calculation is very efficient. 

In order to ascertain the accuracy of the FE-DVR method, the solutions of Eq. (jH]) are 
compared with the solutions obtained by the spectral integral equation method (S-IEM) 



22|, whose accuracy is 1 : 10 as described in Appendix A. The numerical FE-DVR 



solutions are first normalized by comparison with the S-IEM solutions at one chosen radial 
position near the origin, and the error of the normalized FE-DVR function is determined by 
comparison with the S-IEM function at all other radial points r. Since the S-IEM function 
depends on the values of the potential at all points [0 < r < -Rmax], the S-IEM calculation has 
to be carried out to a distance -Rmax large enough so that the contribution from V{r > -Rmax) 
is smaller than the desired accuracy of the S-IEM solution. The same is not the case for 
the FE-DVR solutions ipFE-DVR{f), since the un-normalized solution depends only on the 
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potentials for distances less than r. However, if the normalization of the wave function ([7]) is 
to be accomplished by matching it to sin(fcr) and cos(A;r) at -Rmax in the asymptotic region, 
then the numerical errors that accumulated out to -Rmax will affect the wave function at all 
distances. These errors can be avoided by an iterative procedure for the large distance part 
of the wave function, as will be described in a future publication ^] . 

The results for potentials Vws and Vm are shown figures [3] and HJ respectively. In both 
cases the error of the wave function starts with 10~^^ at the small distances, and increases 
to 10^^° as the distance increases, due to the accumulation of various errors. The accuracy 
for the nonlocal potential K is shown in Fig. [5l The accuracy of the integral (IHD, for a fixed 
size of all partitions as a function of the number of Lobatto points in each partition, is 
shown in Fig. [HI where the open circles represent an upper limit of the estimated accuracy 
as developed in Appendix B, of order (A^ — 2)^. This figure is important because it shows 
the nearly exponential increase of accuracy as A^ increases, until the accumulation of errors 
overwhelms this effect once the value of A^ increases beyond a certain value, 20 for the 
case of Fig. [6l The accuracy of the integral (IHl) for a fixed number A^ per partition, but 
for several different partition sizes, is displayed in Fig. |9l This figure shows that the 
accuracy decreases exponentially with the size of the partition, which can be interpreted as 
an exponential increase of the accuracy with the number of Lobatto points in each partition 
of fixed length. 

Finally, the FE-DVR computing time as a function of the number A^ of Lobatto points in 
each partition is displayed in Fig. [HI where it is also compared with an estimate described in 
Appendix B of the number of floating point operations expected. According to this estimate, 
the time per floating point operation turns out to be ^ 10^^ in a MATLAB computation 
performed on a desktop using an Intel TM2 Quad, with a CPU Q 9950, a frequency of 
2.83 GHz, and a RAM of 8 GB. The dashed line represents the total time required for a 
comparable S-IEM computation. That comparison shows that the FE-DVR method can be 
substantially faster than the S-IEM even though the former has many more support points, 
depending on the radial range and on the accuracy required. Further details are given in 
Table [TTl in Appendix A. 

A comparison between the FE-DVR and a finite difference sixth order Numerov method 
of the accuracy of tan((5) is illustrated in Fig. [71 This comparison shows that for an 
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FIG. 2: (Color on line) The wave functions for the local potentials Vm and VwSi and for the 
nonlocal potential K, described in the text. The wave number is k = 0.5 fm^^ and the potentials 
Vm and Vws are illustrated in Fig. [TJ 
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FIG. 3: (Color on line) Accuracy of the FE-DVR solution of the Schrodinger Eq. for the Woods 
Saxon potential Vws displayed in Fig. [TJ The wave number is A: = 0.5 fm~^, the size of each 
partition is l.fm, and there are 20 Lobatto points per partition. The graph shows the accuracy of 
the wave function ip by displaying the absolute value of the difference between the FE-DVR and 
the S-IEM wave functions. The latter is deemed accurate to 1 : 10^^^. 

accuracy of tan(5) of ^ 10^^, the FE-DVR method requires 15 times fewer meshpoints, and 
is approximately 100 times faster than the Numerov method. More details are presented in 
Appendix C. 
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FIG. 4: (Color on line) The accuracy of the FE-DVR wave function for the potential Vm as obtained 
by comparison with the S-IEM result. The latter is accurate to 1 : 10~^^. The wave number is 
k = 0.5 fm~^, the number of Lobatto points per partition is 20, the size of each partition is 1 fm. 
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FIG. 5: (Color on line) Same as Fig. [3] for the nonlocal Perey-Buck potential K{r,rf). The wave 
number is A; = 0.5/m^^, only one partition was used in the full radial interval from to 15 fm 
using a total of 130 Lobatto grid points. The accuracy of 10^^ is consistent wth the estimate made 
in Eq. (j33p in Appendix B. 

IV. SUMMARY AND CONCLUSIONS 

The accuracy of a hybrid finite element method (FE-DVR) has been examined for the 
solution of the one dimensional Schrodinger equation with scattering boundary conditions. 
This method [l^ uses as basis functions the discrete variable representation Lagrange poly- 
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FIG. 6: (Color on line) Accuracy of the integral J^^^ sin{kr) VM{r) 'ifj{r) dr, obtained with the 
FE-DVR method as a function of the number of Lobatto points in each partition. The length of 
each partition is 1.0 fm, the number of partitions is 100. The potential is Vm, the wave number 
is k = 0.5 fm~^, the accuracy is obtained by comparisom with the S-IEM result which is accurate 
to 1 : 10~". The open circles represent an estimate of the upper bound for the accumulation of 
roundoff errors, given by Eq. (j30p in Appendix B. 
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FIG. 7: (Color on line) This accuracy comparison for tan((5) is performed for the potential Vm and 
k = 0.5 fm^^ in the radial interval [0, 100/m]. The partition sizes in the FE-DVR method have 
a length of 1 fm each, and the number of Lobatto points in each partition is given by 1/100 th 
of the total number of points. Numerov is a 6th order finite difference method with equidistant 
points, as described in Appendix C. 
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FIG. 8: (Color on line) The computing time in MATLAB for the calculations described in Fig. [6l 
The estimate is given by Eq. (j30p . with the factor 10^^^ replaced by 2 x 10~^. The latter represents 
the time for each floating point operation. The dashed line represents the computing time for the 
STEM calculation, described in Fig. [6l 

nomials ii{r), i = 0,1, 2, ..N, on a mesh of Lobatto support points. The accuracy of the 
FE-DVR method is obtained by comparison with a spectral method S-IEM, whose accuracy 
is of the order of 1 : 10^^^. An important advantage of a discrete variable representation 
basis is the ease and accuracy with which integrals can be performed using a Gauss-Lobatto 
integration algorithm that furthermore render the matrix elements {ii{V — E)£j) diagonal. 
This feature permits one to easily solve the Schrodinger Eq. also in the presence of nonlocal 
potentials with a kernel of the form K{r,r'), as is demonstrated in one of our numerical 
examples. Another advantage is that the Galerkin matrix elements of the kinetic energy 
operator T need not be recalculated anew for each partition because they are the same in 
all partitions to within a normalization factor that only depends on the size of the partition. 
A further advantage is that the convergence of the expansion ([9]) with the number N of 
basis functions is exponential, in agreement of what it is the case for bound state finite 
element calculations with Lobatto discretizations [12]. A possible disadvantage may be that 
if the number of the Lagrange polynomials in each partition is very large and/or the number 
of partitions is large, as is the case for long ranged potentials, then the accumulation of 
roundoff and algorithm errors may become unacceptably large. 

In summary, for scattering solutions of the Schrodinger equation the accuracy of the 
FE-DVR method increases exponentially with the number of Lagrange polynomials in each 
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FIG. 9: (Color on line) Accuracy of the integral Jq s'm(kr) VM{r) tp{r) dr, Eq. ([8]), as a function 
of the length of each partition, into which the radial interval [0, 100 fm] is divided. The total 
number = 20 of Lobatto points in each partition is kept constant. The conditions are the same 
as in Fig. [6l This figure shows that the accuracy decreases exponentially with the size of the 
partition. For 20 partitions the computation time is 0.060 s, for 100 partitions it is 0.075 s. 



partition until the accumulation of roundoff and truncation errors overwhelm the result. 
The FE-DVR can easily achieve an accuracy of the order of 10~^° for the scattering phase 
shifts for either local or nonlocal short ranged potentials; it is less complex than the spectral 
S-IEM method but is comparable in the amount of computing time; and, in addition, it is 
substantially more efficient than a finite difference Numerov method. The latter result is 
demonstrated by the fact that the FE-DVR was found to be a hundred times faster than 
the Numerov for an accuracy of 10~® of the scattering phase shift. 

Acknowledgements: One of the authors (GR) is grateful to Professor McCurdy for a 
stimulating conversation on the use of Lagrange polynomials in finite element calculations. 

Appendix A: the S-IEM method ^ 

A version of the spectral method employed here was developed recently [22]. It consists 
in dividing the radial interval into partitions of variable size, and obtaining two independent 
solutions of the Schrodinger Eq. (E]) in each partition , denoted as Y{x) and Z{x). These 
solutions are obtained by transforming Eq. (E]) into an equivalent Lippmann-Schwinger 
integral equation (L-S) and solving the latter by expanding the solution into Chebyshev 
functions, mapped to the interval [—1, +1]. The corresponding discretized matrices are not 
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FIG. 10: (Color on line) The partition distrtibution for the S-IEM method in the radial interval 
[0, 100 fm]. for two different numbers N of the Chebyshev expansion functions in each partition. 
The end point 62 of each partition is shown on the vertical axis, and the corresponding partition 
number is shown on the horizontal axis. The potential is Vm , and the wave number is A; = 0.5 fm^^. 
The accuracy parameter tol in each partition is 10"^^. The computation time for each case is 
approximately the same, 0.2 s,and the accuracy of the wave function in both cases is the same. 
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sparse, but are of small dimension equal to the number of Chebyshev points per partition. 
The solution ip in each partition is obtained by a linear combination of the two independent 
functions Y{x) and Z{x), with coefficients that are determined from the solution of a matrix 
equation of dimension twice as large as the number of partitions, but the corresponding 
matrix is sparse. Details are given in Ref. [22|, and a pedagogical version is found in Ref. 



28|. 

One of the features of the S-IEM method is that the size of each partition is adaptively 
determined such that the accuracy of the functions Y{x) and Z{x) is equal or better than 
a pre- determined accuracy parameter tol, which in the present case is tol = 10^^^. In the 
region where the potential V is small the corresponding partition size is large. When the 
number of Chebyshev expansion functions per partition is large the size of the partitions 
is correspondingly large. As is illustrated in Fig. [10] when is increases from 17 to 33 
the number of partitions decreases from 29 to 6, yet the accuracy of the respective wave 
functions is approximately the same, 1 : 10~^^, and the computing time is also approximately 
the same, 0.2s. 
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FIG. 11: (Color on line) The y-axis illustrates the absolute value of the difference between two 
S-IEM wave functions, calculated with accuracy parameters tol = 10^^^ and 10~^^, respectively, 
for potential Vm and k = 0.5 fm~^. This difference is less than 4 x 10^^^ for all values of r. 

For the present S-IEM benchmark calculations the value of is 17, and for the case of 
Vm the maximum value of r is 100 fm. Such a large value is required because the potential 
decays slowly with distance and becomes less in magnitude than 5 x 10~^^ only beyond 
r = 100 fm. Had the potential been truncated at a smaller value of r, then the truncation 
error would have propagated into all values of the wave function and rendered it less accurate. 
The accuracy of the S-IEM wave function can be seen from Fig. [TTl which compares two 
S-IEM wave functions with accuracy parameters tol = 10^^^ and 10^^^, respectively. The 
result is that the accuracy of the lEM wave function for = 17 and -Rmax = 100/m and 
tol = 10"^^ is 4 X 10~^\ and that for tol = 10~^^ the accuracy is better than 10^^^. 

The wave functions are normalized such that their asymptotic value is given by Eq. ([7]). 
The corresponding values of tan((5), Eq. ([H]), for potentials Vm and Vws and a wave number 
k = 0.5 /m-i are 2.6994702502 and -1.7107344227, respectively. Table [H shows the number 
of partitions, the accuracy of tan((5), and the computing time of the S-IEM method for 
various tolerance parameters inputted into the code for the potential Vm, with k = 0.5. 
The number of Chebyshev polynomials in each partition is 17, the total number of points 
displayed in the third column is equal to 17 times the number of partitions. The error of 
tan(5) is obtained by comparing the value of tan(5) for a particular tolerance parameter 
with the value obtained for tol = 10^^^. 
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Tol. 


Part'ns 


Points 


Err[taii{5)] 


time {s) 


10-12 


37 


629 




0.178 


10-10 


25 


425 


4.6 X 10-12 


0.181 


10-8 


17 


289 


7.7 X 10-11 


0.171 


10-6 


11 


187 


5.2 X lO-'^ 


0.165 


10-4 


7 


119 


2.8 X 10-4 


0.162 


10-2 


5 


85 


6.5 X 10-2 


0.161 



TABLE I: Accuracy and computing time for the S-IEM method 



of Pts. 


err [tan((5)] 


time{s) 


2000 


10-10 


0.075 


1300 


10-8 


0.050 


1200 


10-6 


0.047 


1000 


10-4 


0.045 


700 


10-2 


0.042 



TABLE IL Accuracy and computing time for the FEM-DVR method 



For the case of a nonlocal potential K the division of the radial interval into partitions 
is not made because the effect of the nonlocal potential would extend into more than one 
partition, making the programming more cumbersome. For the case of a kernel K{r,r'), 
described in Ref. [25|, the accuracy of the S-IEM result [25,] is also good to 1 : IQ-n, as is 



shown in Fig. 7 of Ref. {25 1. 

For comparison with the S-IEM method some characteristics of the FE-DVR method are 
shown in Table [Tll The potential and the wave number is the same as in Table HI the radial 
interval [0, 100 fm] is divided into 100 partitions of length of 1 fm each, and the number 
of Lobatto points per partition in all partitions is the same but is progressively varied from 
20 to 7, as shown in the table. If one compares the entries in table H] with those in |Tl] that 
correspond to approximately the same accuracy of 10^ for tan(5) one notices that the FE- 
DVR method needs approximately 7 times more support points than the S-IEM, yet the 
computing time is between 2 and 3 times less. This remark attests to the efficiency of the 
FE-DVR method. 
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Appendix B: The round off errors in the FE- 
DVR method. 

The notation is as follows: is the number of Lobatto points in each partition, which is 
also equal to the number of Lagrange polynomials in each partition, and Np is the number 
of partitions. The largest contribution to the roundoff errors is expected to arise from the 
solution of Eq. fl24l) for the N — 2 expansion coefficients. This equation is of the type 
M/3 = 6, where /3 is the column vector of the N — 2 expansion coefficients and M is a 
matrix of dimension N — 2, whose solution requires 4 x (A^ — 2)^ floating point operations. 
For the case that the floating point roundoff error of the computer is e and that the errors 
accumulate linearly, an upper bound for the total error Et is 

eT^^*Np{N -2)^e. (30) 

For Np = 100 and for e = 10~^^, which is the value for the calculations done in MATLAB, 
one obtains an upper bound for the values of Et that are plotted in Fig. E] as a function of 

N 

~ 4 X 10-^^(A^ - 2fe. (31) 

The floating point error that occurs in the calculation of the Lagrange functions ii{x) 
is much smaller. The numerator contains A^ factors x — Xi , each of which can be written 
as Aj + e, where Aj is proportional to the length of each partition. Hence the error of the 
product is ~ A^ + NA^~'^e, where A is an average value of the x — Xi. A similar argument 
holds for the denominator and if the error of the numerator adds linearly to the error of the 
denominator, then an upper bound for the total error of a Lagrange function is ~ 2Ne/A. 
This is much less than the error in Eq. (!30l) . 

For the case of the nonlocal calculation the conditions above are different. There 
is only one partition of length L = 15 fm, the number of Lobatto points is 
130, and the order of each polynomial i{r) is 129. The error in the calculation 
of the Lagrange polynomials, or their derivatives at each meshpoint, is < 2(A^ — 
l)e/A. Since the error in the calculation of the matrix element of {(P/dr"^) has 
A^ terms according to Eq. ([5]) that could lead to an error of ~ 2N{N — 
1)6/ A, assuming that all e errors add linearly. The solution of of Eq. (12^ requires 4(A^ — 
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2)^ operations and thus an upper bound of the total hnear accumulation of e errors is 

~ 2N{N - l){e/L)A{N - 2f = 1.8 x 10~^ (32) 

Since the e errors do not accumulate linearly, the expected upper bound for the error could 
be 

< 2{N - l){e/L)A{N - 2)^ = 1.4 x 10"^ (33) 
The above estimate is consistent with the accuracy found numerically in Fig. |5l 

Appendix C: Comparison with a finite difference 
method. 

The finite difference method used for this comparison is Milne's corrector method, also 
denoted as the Numerov method, given by Eq. 25.5.21 in Ref. 23|. In this method the 
error of the propagation of the wave function from two previous points to the next point is 
of order h^, where h is the radial distance between the consecutive equispaced points. The 
calculation is done for the potential Vm and for k = 0.5 fm~^ as follows. 

A value of h is selected and the Milne wave function is calculated starting at the two initial 
points r = h and r = 2/i by a power series expansion of the wave function for the potential 
Vm- The values of the wave function for the additional points 3h,Ah, ...are obtained from 
Milne's method out to the point r = 20/m. The wave function is normalized to the S-IEM 
value at r = 2 fm, and the error at r = 18 fm is obtained by comparison with the S-IEM 
value at that point. The result for a sequence of h values is illustrated in Fig. [121 For each 
value of h the wave function is calculated out to r = 100 fm by Numerov's method, and 
;he integral ([8]) is calculated by the extended Simpson's rule, given by Eq. (25.4.6) in Ref. 



23| . The error is determined by comparison with the S-IEM result 2.6994702502 for tan(5). 
A comparison with the FE-DVR method is shown in Fig. [7] and Table IIIII displays the ratio 
Numerov/FE-DVR of the total number of points and of the time of the two methods for two 
accuracies of tan(5). . More detail of the error and the computing time for the Numerov 
method is displayed in Table HVl 

. The calculation is done in MATLAB performed on a desktop using an Intel TM2 Quad, 
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10 10 10 10 

Total number of points 



FIG. 12: (Color on line) The error of the Numerov wave function at r = 18 fm, as a function of 
the number N of meshpoints in the interval [0,20/m]. The distance h between points is 20 /N. 
For each h the wave function is normalized to the S-IEM wave function at r = 2 fm. The wave 
number is k = 0.5 fm^^, the potential is Vm- 



accuracy tan(5) 


ratio of N of Pts 


time ratio 


10-6 


~ 1 


20 


io-« 


15 


100 



TABLE III: The Numerov/FE-DVR ratio of the required total number of meshpoints and the 
respective computational times 



with a CPU Q 9950, a frequency of 2.83 GHz, and a RAM of 8 GB. 
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